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1.   Introduction.   Multistep  methods  are  the  oldest  and  the  most  important 

class  of  numerical  methods  for  solving  systems  of  ordinary  differential  equations, 

Implementations  of  these  methods  have  become  increasingly  sophisticated 

in  the  last  two  decades.   One  paper  having  a  considerable  impact  on  the 

development  of  these  algorithms  is  that  of  Nordsieck  [1962],  which 

introduces  a  class  of  formulas  closely  related  to  multistep  formulas. 

It  is  the  purpose  of  the  present  paper  to  explore  thoroughly  the 

relationship  between  multistep  and  Nordsieck  formulas  on  a  uniform  mesh. 

The  results  apply  to  a  limited  extent  to  variable  order  variable  stepsize 

algorithms  which  discourage  order  and  stepsize  changes,  because  such 

algorithms  often  produce  sequences  of  values  computed  with  the  same 

formula  and  stepsize. 

Consider  the  system  y'  =  f(t,  y)  with  exact  solution  y(t). 

Given  starting  values  y.,  y! ,  j  =  0(l)k-l,  a  linear  k-step  formula  on  a 

mesh  t  :=  t_  +  nh,  n  =  1(1)N,  determines  approximations  {y  }  to  the  values 
n    0  n 

{y(t  )}  by  means  of  p(E)y  .  =  ha(E)y'  .  and  y'  =  f(t  ,  y  )  where 
n  n-k         n-k      n      n   n 

k  k 

p(£)  :=  Z     a.   C   J  ,  a(£):=  £  6.  ?  3  , 
j=0  J  j=0  J 

2    2 
and  E  is  the  shift  operator.   It  is  assumed  that  ex  ^  0  and  a,  +3,  >  0. 

0  k    k 

For  convenience  a  normalization  such  as  a  =  1  or  Z$  =  1  is  not  imposed. 
If  3n  $   0,  the  formula  is  implicit  and  requires  the  solution  of  a  system 
of  nonlinear  equations,  which  is  always  possible  for  small  enough  h  if 
f(t,  y)  is  Lipschitz  continuous  as  a  function  of  y.   In  §2  a  simple 
reformulation  of  the  general  k-step  formula  is  given  which  requires  that 
only  k  values  be  saved  between  steps.   In  practice,  k+m  values  might  be 
stored,  from  which  it  is  possible  to  construct  a  predictor  formula  of 
order  k+m-1  except  for  pathological  cases  where  p(£)  and  a(£)  have  a 
common  factor. 


In  §3  linear  multistep  formulas  are  described  in  terms  of  a 

local  polynomial  approximation  p  (t)  which  interpolates  the  k-hn  values 

which  are  retained  after  completion  of  the  n-th  step.   It  is  from  this 

basis  free  description  that  equivalent  forms  of  multistep  formulas  are 

derived.   One  convenient  representation  of  the  polynomial  p  (t)  is  in 

n 

terms  of  the  scaled  derivatives  h  p  J  (t  )/j!,  j  =  0(l)k+m-l.   Multistep 

formulas  represented  in  this  way  belong  to  the  class  of  formulas 

introduced  by  Nordsieck.   Another  is  in  terms  of  ordinates  p  (t   .), 

n  n-j 

i  =  0(l)k  -1,  and  derivatives  hp'(t   .),  j  =  0(l)ko-l,  where  kn  +  k.  = 
1  n  n-j  2  12 

k+m,  which  yields  the  modified  multistep  formulas  of  Gear  [1971a,  p.  150] 

Backward  differences  of  ordinates  or  derivatives  is  also  possible. 

General  linear  Nordsieck  formulas  compute  approximations 

a  =:  [y  ,  hy',...,  h  y    /q!l   to  the  vector  of  scaled  derivatives 
~n     n    n         n 

[y(t  ),  hy'(t  ),...,  hqy^(t  )/q!]   by  the  recurrence 
n        n  n 

a  :=  Pa  .  +  £6 
~n    ~n-l   ~  n 

where   6      is   chosen  so   that  y'    =   f(t    ,   y   );    that   is,    6      is   implicitly 
n  n  n       n  n 


defined  by 


e.Pa     ,    +  JL6     =  hf(t    ,    e_Pa     .    +  ln6   )    . 
~1   ~n-l  In  n     ~0  ~n-l  0  n 

Here  P   is   the  Pascal   triangle  matrix  defined  by 

P.     :-    (.h,  0<i<j<q, 

T  T 

e.  is  the  i-th  unit  vector,  and  I   =:  [£.,  £,,...,£]   is  a  vector  of 
~j  0   1       q 

parameters  which  characterize  the  formula.   It  is  assumed  that  £  ^  0. 
Note  that  zero-origin  indexing  is  being  used.   Nordsieck  [1962]  derives 
particular  choices  for  %   which  optimize  stability  and  accuracy  for 
q=2,  3,  4,  5,  6,  and  he  implements  the  formula  for  q  =  5  as  a  variable- 
step  computer  program  for  the  Illiac  I.   By  direct  calculation  he  shows 


that  there  is  a  connection  between  this  formula  and  the  sixth  order 

Adams-Moulton  formula;  specifically,  it  is  shown  that  the  zero-th 

and  first  component  of  a  ,  a   ,,...,  a   c  satisfy  the  difference  equation 
r  ~n  ~n-l      ~n-5 

of  the  Adams-Moulton  formula. 

Further  understanding  of  Nordsieck  formulas  resulted  from 

the  work  of  Descloux  [1963].   An  actual  equivalence  is  shown  between 

Nordsieck' s  (q+l)-value  formula  and  the  (q+l)-st  order  Adams-Moulton 

formula,  in  the  sense  that  a  linear  transformation  was  constructed 

relating  the  Nordsieck  vector  to  the  vector  of  values  [y  ,  hy',..., 

n    n 

T 
hy'   , , ]  .   An  equivalence  was  also  shown  for  the  q-th  order  Adams-Moulton 
n-q+1 

formula  and  a  (q+l)-value  Nordsieck  formula  with  a  slightly  different 
£  vector.   This  second  family  of  formulas  is  used  in  the  nonstiff  option 

of  DIFSUB  (Gear  [1971b])  and  its  descendants.   Furthermore,  Descloux 

T 
showed  that  in  both  cases  the  predicted  value  y  _:=  e,_Pa   ,  satisfies 

n,U   ~U  ~n-l 

the  difference  equation  of  the  k-th  order  Adams-Bashf orth  formula. 

General  k-step  formulas  were  also  considered  and  equivalent  Nordsieck 

formulas  were  constructed  having  the  more  general  form 

(1.1)  a  =  Aa   _  +  £6 

~n    ~n-l   ~  n 

where  the  dimensionality  of  the  arrays  might  be  as  high  as  2k  and  where  A 

is  not  necessarily  the  Pascal  triangle  matrix.   In  the  appendix  of  Descloux' s 

report  is  a  theorem  stating  that  the  values  y  and  f  computed  by  formulas 

n      n 

even  more  general  than  (1.1)  satisfy  the  difference  equation  of  an 
associated  linear  multistep  formula.   The  proof  is  not  given  by  Descloux 
because  it  "is  rather  painful  and  long."  However,  in  a  later  paper  of 
Osborne  [1966]  a  similar  result  is  stated  and  a  short  and  elegant  proof 
is  given. 


In  §4  a  one-to-one  correspondence  is  established  between  (i) 
the  class  of  all  (q+l)-value  linear  Nordsieck  formulas  and  (ii)  the 
class  of  all  linear  multistep  formulas  of  formal  order  at  least  q 
and  stepnumber  at  most  q.   This  correspondence  is  given  in  closed  form 
by 

^"VO  -  det(Sl-P)  Jai-P)"1*, 

and 

Cm_1a(5)  =  det(£l-P)  eJcCl-P)"1^  . 

Osborne  [1966]  gives  an  expression  for  E,        p(£)/(£-l)  similar  to  that 

given  above  and  proves  a  correspondence  between  Nordsieck  and  multistep 

formulas  of  order  at  least  q+1  (although  it  is  not  clear  that  I     is 

correctly  chosen  in  the  proof).   Presumably  at  that  time  only  formulas 

of  optimal  order  were  considered  to  be  of  practical  value.   More 

recently  Wallace  and  Gupta  [1973]  obtain  expressions  for  E,        p(£)  and 

E,       o(E,)   which  are  equivalent  to  those  given  above. 

In  §5  it  is  shown  that  the  zero-th  and  first  components  of 

a  ,  a   ..,...,  a  .  satisfy  the  difference  equation  of  the  corresponding 
~n  ~n-l      ~n-k 

multistep  formulas.   Also,  if  p(£)  and  o(E,)   have  no  common  factor,  then 

the  Nordsieck  formula  is  shown  to  be  equivalent  in  the  sense  of  Gear 

[1971a,  p.  143]  to  the  corresponding  multistep  formula,  which  means  that 

a   is  a  linear  transformation  of  certain  linear  combinations  of  the 
~n 

values  y   . ,  y1  .,  i  =  0(l)k-l.   It  is  proved  that  there  is  no  equivalence 
yn-j   •'n-j 

to  a  multistep  formula  if  p(?)  and  a (£)  have  a  common  factor. 

In  §6  predictor-corrector  Nordsieck  formulas  are  considered. 
An  equivalence  to  predictor-corrector  multistep  formulas  is  shown  to 
be  the  case  for  P(EC)*  formulas  but  not  for  P(EC)*E  formulas. 

In  §7  recent  applications  of  equivalence  results  are  discussed. 


2.   Minimum  Storage  for  Multistep  Formulas.   The  computation  of  each  solution 

value  y  and  its  derivative  y'  by  a  general  linear  multistep  formulas  would 
^n  "n 

seem  to  require  saving  2k  previously  computed  values  y   ,,y   „,...,  y  . 

n-i   n-z       n-k 

and  y'  ..  ,  y'  _,...,  y'  .  .   For  this  reason  it  has  been  considered  desirable 
n— i   n— 2       n-k 

to  have  formulas  with  most  of  either  the  a  or  3  coefficients  set  to  zero 
(see,  for  example,  the  "minimum  storage  methods"  of  Dill  and  Gear  [1971]). 
However,  there  is  an  easy  way  to  get  by  with  just  k  values  without  re- 
evaluating the  function  f.   Define 

sJ  :=  I      (-a.  ...y   .  +  hR   ...y'  .) 
n    .  n  k-i+1  n-i     k-j+1  n-i 

i=0      J 

for  j  =  0(l)k-l.   In  terms  of  the  partial  sums  s  ,  the  linear  multistep 
formula  becomes 

(2.1)      solve  ay  =  hg  f(t  ,  y  )  +  s  "   for  y   , 
U  n     (J   n   n     n-i      n 

set     y'  =  f(t  ,  y  )  , 
n      n   n 

Sr,      =   "aiy„  +  h3l^  +  S„      1   ' 

n       In     In    n-1 

s   =  -a   ny  +  h3   y1  +  s     , 
n     k-1  n     k-1  n    n-1 

s  =  -a.  y  +  h8  y '  . 
n     k  n     k  n 

Thus  only  the  k  partial  sums  s   -,,s   .,,...,  s   ,  need  to  be  saved  in  order 

n-1   n-1       n-1 

to  advance  the  computation.   Of  course,  this  is  not  the  only  set  of  values 

that  could  be  saved;  any  nonsingular  linear  transformation  of  them  could 

be  used  instead. 

In  practice  the  implicit  equation  (2.1)  is  approximately  solved 

by  some  iterative  process.   An  initial  guess  y    is  obtained  with  a 

n,  0 

predictor  formula 


k*  k* 

a*y    +  I     a*y   .  =  h  Z  3*y' 
0  n,0        j'n-j  j'n-j 

and  subsequent  iterates  y    are  given  by  an  iteration  of  the  form 

n,u 

y    ,  :=  y    -[anl  -  h3„J    ]_1  [a„y    -  hgrtf    -  sk"^]  , 
^n,y+l   'n,y    0      0  n,y     l  (Tn,y     0  n,y    n-1   ' 

y  =  0(l)M(n)-l  , 

where  f    :=  f(t  ,y   ) ,  1  is  the  identity  matrix  and  J    is  a  matrix 
n,y       n  n,y  n,y 

depending  on  the  type  of  corrector  iteration.   For  functional  iteration 

J    :=  0,  for  the  Newton  iteration  J    :=  f  (t  y   ),  and  for  the 
n,y  n,y     y  n/n,y" 

chord  iteration  J     :=  f  (t  ,y   .,)  .   The  number  of  iterations  M(n) 
n,y     y  n  ^n,0  J 

might  vary  from  step  to  step.   After  determining  y   :=  y  w,  .  a  final 

n     n,M(n) 

function  evaluation  may  or  may  not  be  performed.   For  a  P(EC)*  formula 

y'  :=  f  ...   N  ,  is  accepted  and  for  a  P(EC)*E  formula  y'  :=  f  =  f  w/  .. 
^n     n,M(n)-l        r  Jn  n    n,M(n) 

For  our  purposes  the  essential  difference  between  the  two  types  of 

predictor-corrector  schemes  is  that  the  values  y  and  y'  satisfy  the 

n      n 

corrector  formula  for  the  P(EC)*  scheme  and  they  satisfy  the  differential 

equation  for  the  P(EC)*E  scheme.   In  fact,  one  could  view  the  P(EC)* 

formula  as  computing  values  y  ,  y'  that  satisfy  a  slightly  different 

n   n 

differential  equation,  which  implies  that  equivalence  results  for  linear 
formulas  apply  equally  well  to  P(EC)*  formulas. 

How  many  values  must  be  saved  between  steps  for  these  predictor- 
corrector  formulas?   The  answer  would  seem  to  be  the  row  rank  of  the 
following  matrix: 


aT 


a 


K-l 


a. 


a* 
K 


a* 
K-l 


a. 


a. 


°K 


a 


K 


K 


K-l 


PK 


6* 
PK-1 


K 


ft* 
PK 


K 


a* 

a* 

•     "a*. 

3* 

RA 

•     3* 

1 

2 

K 

1 

2      * 

K 

where  K  :=  max{k,  k*}.   The  row  rank  of  this  matrix  is  at  least  K  and 
could  be  as  high  as  k+k*.   For  reasons  of  storage  economy  the  predictor 
should  be  chosen  so  that  the  rank  of  the  matrix  is  low.   The  most 
economical  choice  would  be  to  use  a  predictor  of  the  form 

0        1       L  J       k-l 

V0  =  Vn-l  +  Vn-l+-+Tk-lSn-l 

where  the  coefficients  Y.  are  chosen  to  make  the  formula  of  order  k-l 
or  greater. 

A  predictor  of  order  k-l  may  not  be  accurate  enough  for  a 
corrector  of  order  k  or  greater.   For  this  reason  it  may  be  preferable 
to  store  more  than  k  values  between  steps.   For  a  predictor  of  order  at 
least  q  where  q  >_  k-l,  one  would  expect  to  need  an  additional  m  =  q-(k-l) 
values.   This  can  be  done  by  including  information  from  preceding  meshpoints 
Since  each  step  of  the  computation  introduces  exactly  one  new  item  of 
information,  namely  y1 ,  the  following  set  of  k+m  values  is  suggested: 


n-m 


j  =  0(l)k-l  , 


hy*  .  ,   j  =  0(l)m-l  . 


From  these  values  one  can  apply  the  corrector  formula  to  generate 


y    .. ,  8    ,.,...,  s    ,.»  j  "  l(l)m.   Moreover,  if  the  corrector 
n-m+j    n-mf  j       n-m+j 

formula  is  of  order  q  or  more,  then  the  polynomial  of  degree  q  or  less 
which  interpolates  the  original  set  of  k+m  values  also  interpolates 
the  values  generated  from  them.   An  alternative  set  of  values  from 
which  the  others  can  be  recovered  (by  reversing  the  forward  step 
procedure)  is  the  following: 

sj   ,   j  =  0(l)k-m-l*  , 

n-m 

yn-j'   j  =  °*(1)m-1  ' 

hy'  .,   j  =  0(l)m-l  , 
n-j 

where  0*  -  0,  1*  =  1  if  m  <   k  and  0*  =  1,  1*  =  0  if  m  =  k+1.   This  set 

has  the  important  advantage  of  requiring  values  from  only  k+0*  meshpoints 

instead  of  k+m  meshpoints,  although  it  is  a  little  more  awkward  to  specify, 

The  case  m  =  k+1  can  only  apply  to  the  trapezoid  formula,  the  Milne 

formula,  and  other  (unstable)  k-step  formulas  of  maximal  order  2k.   For 

this  case  the  set  of  values  y   .,  i  =  l(l)k,and  hy'  .,  i  =  0(l)k, is  used 

n-j   J  n-j 

to  define  a  polynomial  of  degree  at  most  2k.   This  polynomial  interpolates 

y  because  y  is  determined  from  the  other  values  by  a  corrector  formula 
n         n 

which  is  exact  for  polynomials  of  degree  at  most  2k. 

As  an  example,  a  k-th  order  predictor  for  the  (k+l)-st  order 
Adams-Moulton  formula  would  be  constructed  from  the  values 

j 

.  _  k-i+i  n-l-i  n        n 

i=0    J 

or  equivalently  linear  combinations  of  these  values:   y  ,  and  hy'  ., 

n       Jn-2 

j  =  0(l)k-l.   The  predictor  so  constructed  is  the  k-th  order  Adams- 
Bashforth  formula. 


It  is  not  obvious,  however,  that  such  a  predictor  can  always 

be  constructed.   Clearly  the  question  is  of  practical  interest  only  for 

implicit  formulas  (p,  a)  of  order  at  least  q.   It  is  claimed  that  the 

construction  of  a  suitable  predictor  is  possible  if  and  only  if  a  unique 

polynomial  p(t)  of  degree  q  or  less  is  uniquely  determined  from  the 

values  L^p(t   ),  j  =  0(l)k-l,  and  p'(t   . ),  j  =  0(l)m-l,  where  the 
h   n-m  n-j 

operators  L?    are  defined  by 

i        j 
L^z(t)  =  Z   (-a.  .^z(t-ih)  +  h6,  .,.z'(t-ih))  . 
i=0    k_J+:L  k-j+i 

Clearly,  the  possibility  of  constructing  p(t)  implies  the  existence  of 

a  suitable  predictor.   On  the  other  hand,  if  p(t  ,,)  can  be  determined 

n+1 

by  a  predictor,  then  from  the  corrector  formula  we  can  get  p'(t   . ) . 

This  can  be  repeated  to  generate  p(t  l0),  p(t  ,_),...  until  there  are 

n+z      n+3 

enough  values  to  construct  p(t)  by  interpolation  of  ordinate  values. 
Necessary  and  sufficient  conditions  for  the  possible  construction  of 
p(t)  is  given  by  the  theorem  that  follows.   First  a  precise  definition 
of  order  is  needed. 

Definition.      A  linear  multistep  formula  (p,  a)  is  of  order 
at  least  q  if 

Mi±f}  .  iog(i+z)  +  o(zs+1) 

and  is  of  formal  order   at  least  q  if 

p(l+z)  =  log(l+z)  a(l+z)  +  0(zq+1) 

(cf.  Gear  [1971a,  p.  119]). 

THEOREM  2.1.     Let   (p,   a)  be  a  linear  k-step  formula  of  order 

at   least  q.      Values  L^p(t        ),   j   =  0(l)k-l,   and  p'(t      .),   j   =  0(l)m-l, 

h   n-m  n-j 

uniquely  determine  a  polynomial   p(t)  of  degree  at  most   q  if  and  only 
if  p(£)  and   a(^)  have  no  common  factors. 


10 

Proof.      If  polynomials  are  identified  with  the  column  vector 
of  their  scaled  derivatives,  then  the  unique  existence  of  p(t)  is 

equivalent  to  the  linear  independence  of  the  rows  of  the  matrix 

i  nT   -kT      ^T     T   T  -1       Tnl-nK   . 
T  :=  col{A0,  A  ,...,  ^    »  e  ,  e  P   ,...,  e^V        )  where 

,T     ^   r        T  ,  0      Ti_— i— m 
A.  :=  Z  {-a     e  +  3,  ...e.JP 
~J    ±=0    k-j+i~0    k-j+i~l 

and  P  is  the  (q+1)  *  (q+1)  Pascal  triangle  matrix.   The  row  vectors 

T 
A  ,  j  =  0(l)k-l,  have  a  simple  generating  function,  for 

(-p(S)eJ  +  a(?)e^)(I-?P"1)"Vm 

=  Z   Ck"j(-a.e^  +  B.e*)   Z  6"1"1 
j=0        J~U    J    i=0 

k"1   £  T 
£=0 

T  k-£    T  - 

where  we  have  used  the  fact  that  A  P    =  0  .   Let  p(£)  =  tt(£)  p(£)  and 

a(£)  =  tt(D  a(0  where  tt(^)  is  the  g.c.d.  of  p(£)  and  a(^).   Then 

—   —  —     T—     T— 1— 1  — m 

(p,  a)  is  of  order  q  or  greater,  and  so  (-  p(^)e  +  o(E,)e   )  (1-^,?     )   P 

-T 
is  a  generating  function  for  the  y.»  which  are  defined  in  the  obvious 

way.   The  relation 

(-P(C)eJ  +  cr(DeJ)(I-£P~1rV",n 

=  Tr(0(-P(0ej  +  a(?)eJ)(I-?P"1)"1P-m 

implies  that  the  first  k  rows  of  T  are  linear  combinations  of  the  row 

-T 
vectors  y.»  J  =  0(l)k-l.   If  p(£)  and  a(£)  have  a  common  factor  so  that 

k  <  k,  then  p(t)  is  obviously  not  uniquely  determined.   Assume  that 

p(£)  and  a(£)  have  no  common  factors.   Suppose  that  the  rows  of  T  are 

linearly  dependent.   Then  Tv  =  0  for  some  v  ^  0.   Hence 
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(-p(Oel  +   c(Oe*)(I-£P  X)  XP  "v  =  0 

or  equivalently 

(2.2)  P(C)tq(?)  =  a(?)T1(C) 

where  T  (£)  =  det ( (I-?P~1)Pin)eT(I-?P~1)~  P~mv,  which  by  Cramer's  rule  is 
simply  the  determinant  of  (I-£P  )P  with  the  i-th  column  replaced  by  v. 
Thus  T . (O  is  a  nonzero  polynomial  of  degree  q  or  less.   Furthermore 

^(o  =  -a-^+V^u-rVV^v 

oo 

=  (1-C)     £  5    e ,P     v 

j-o 


=  (-l)^C   e  Pv  +  lower  order  terms 


and  so  t.. (?)  is  of  degree  at  most  k-1.   Hence  from  (2.2)  at  least  one  of 

the  k  roots  of  p(?)  must  be  a  root  of  G(£).   This  is  a  contradiction  and 

therefore  the  rows  of  T  are  linearly  independent. □ 

COROLLARY.  The  values   L^p(t   ),  j  =  0(l)k-m-l*,  p(t   .), 

n   n-m  n-j 

j  =  0*(l)m-l,  and   p'(t  _.),  j  =  0(l)m-l,  uniquely  determine   p(t)  if  and 
only  if   p(£)  and   a(?)  have  no  common  factor. 

Remark   1.   It  is  not  enough  to  have  formal  order  of  at  least 

2 
q.   Consider  p(£)  =  (£-1)   and  0(O    =  0;  the  values  -p(t  )  and  2p(t  )  - 

n         n 

p(t   ., )  uniquely  determine  a  polynomial  of  degree  one  or  less. 
n-1 

Remark   2.   Theorem  2.1  appears  to  be  related  to  a  result  of 
Dahlquist  [1975]  concerning  the  equivalence  between  a  linear  k-step 
formula  and  the  corresponding  one-leg  k-step  formula  under  the  assumption 
that  p(?)  and  a(£)  have  no  common  factor. 

Apart  from  the  need  for  a  predictor  there  is  another  complication 
that  affects  the  amount  of  storage  needed.   In  practice  it  has  been  found 
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desirable  to  vary  the  stepsize,  and  thus  a  variable  step  form  of  the 
multistep  formula  must  be  used.   For  any  given  fixed  step  formula  there 
are  numerous  variable  step  formulas,  some  of  which  require  less  storage 
than  others.   This  topic  will  be  studied  in  a  future  paper. 
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3.   Construction  of  Linear  Nordsieck  Formulas.   In  this  section  linear 
multistep  formulas  are  described  in  terms  of  local  polynomial  approximations, 
and  an  equivalent  Nordsieck  formula  is  constructed  by  representing  the 
polynomials  in  terms  of  their  coordinates  with  respect  to  local  monomial 
bases. 

Consider  a  linear  k-step  formula  of  order  at  least  q  where 
q  >_  k  such  that  p(£)  and  G (£)  have  no  common  factor.   (The  case  q  =  k-1 
is  treated  at  the  end  of  this  section.)   Then  as  a  consequence  of  the 
corollary  to  Theorem  2.1  a  unique  polynomial  of  degree  at  most  q  is 
determined  by  the  q+1  conditions 

Lk-l(tn-l-J  "  Sn-1-B-     i   -   «»*-l'. 
P„-l(t„-l-j)="yn-l-j'      3-0*(Dm-l, 
K-lK-l-S  =  K-i-y  J-0(D-1. 

The  polynomial  p   , (t)  can  be  regarded  as  an  approximation  to  the  solution 
n-1 

y(t)  near  t  =  t    ,  and  it  can  be  used  to  obtain  a  predicted  value 
n— 1 

y   _  =  p   ,(t  )  as  an  initial  approximation  to  y  .   In  terms  of  the  values 
n ,  0    n-  In  n 

y«  i y„  ir»  and  y'   ,...,  y'   ,  the  relation  y    =  p    (t  )  becomes 

n-1       n-k       n-1       n-k  n,0    n-1  n 

an  explicit  linear  k-step  formula. 

Advancing  the  numerical  solution  one  step  yields  values  which 

determine  the  polynomial  p  (t).   However,  there  is  a  more  direct  way  of 

n 

expressing  p  (t)  in  terms  of  p   ,(t).   First,  we  examine  how  y  and  y1 
n  n-1  n      n 

are  determined.   From  the  multistep  formula  we  have 

ay  =  h3ny'  +  s  ~ 
On     On    n-1 

and  because  L.  annihilates  p   , (t)  we  have 
h  n-1 
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any«  n  =  heny'  n  +  Lu~  P   i  <c   i> 
0  n,0     0  n,0    h   n-1  n-1 

where  y'  _  =  p'  1 (t  ) .   It  follows  from  the  defining  conditions  that 
n,U    n-1  n 

k-1  k-1 

K  p„  i(t«  i^  =  s„  T   Hence  an(yn~yn  r)  =  hen(yJ,_yJ,  J-   Also 

n   n-i  n-l     n-1         U  n  n,U      U  n  n,U 

y'  =  f(t  ,  y  ).   Together  these  last  two  equations  define  y  and  y'. 
n      n   n  n      n 

A  more  convenient  way  of  expressing  this  is  to  introduce  an  increment 

6  which  satisfies 
n 

hy'  n  +  a_6  =  hf (t  ,  y    +  0  6  ) 
n,0    On       n   n,0    0  n 

and  set 

yn  =  yn,0  +  306n  ' 

hy'  =  hy'   +  a_6   . 
n     n,0    0  n 

Second,  we  construct  d  (t)  :=  p  (t)  -  p   ,(t).   For  j  =  0(l)k-m-l* 

n        n       n-1 

Lhdn(tn    J     =    LhPn(tn     J     "     {~\     ^r,     1(tn    J     +    h3V     iK     1  (t«    J 

h  n     n-m  h  n     n-m  k-j   n-1     n-m  k-j   n-1     n-m 

+  Lrj"1?      .(t  )} 

h       n-1     n-l-m 

=  sj       -   {-a.     .y         +  hR      .y'       +  sj"^     } 
n-m  k-j   n-m  k-j    n-m         n-l-m 

=  0    ; 

for  j    =   0*(l)m-l 

jVn   if   3    =   0    , 
dn(tn_.)   =< 

3  1    0        otherwise; 


and   for  j    =  0(l)m-l 


Hence 


jVn   if   J    =   °    ' 

hdn(tn-J    =\ 

I    0       otherwise. 


t-t 

d    (t)   =  A(-r^)6 
n  h        n 
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where  A(x)  is  the  unique  polynomial  of  degree  at  most  q  which  satisfies 

L^A(-m)  -  0  ,    j  -  0(l)k-m-l*  , 

30  ,    j  =  0*(1)0  , 
j  =  l(l)m-l  , 

I  an     •    n 

A'(-j)  -<  °  «  J  =°  ' 

\J)     ,      j  =  l(l)m-l  . 

Briefly,  we  have  shown  that  the  numerical  solution  is  advanced  one 

step  by  setting 

t-t 

(3.1)  p  (t)  =  p    (t)  +  A (-t-^)6 

n       n-l         n    n 

where  6   is  chosen  so  that  p  (t)  satisfies  the  differential  equation  and 
n  n 

A(x)  is  characteristic  of  the  multistep  formula. 

The  polynomial  formulation  of  the  k-th  and  (k+l)-st  order 
Adams-Moulton  formulas  was  discovered  by  Descloux  [1963].   Schemes  based 
on  general  choices  of  A(x)  are  discussed  by  Skeel  [1973]  and  by  Wallace 
and  Gupta  [1973],  who  give  an  interesting  interpretation  of  polynomial 
schemes  in  terms  of  polynomial  predictive  filters.   They  derive  new 
formulas  for  stiff  problems  by  choosing  A(x)  to  be  a  monic  polynomial 
which  best  approximates  zero  for  x  <  0.   Different  types  of 
approximations  yield  different  formulas.   Still  more  formulas  are 
presented  in  Gupta  and  Wallace  [1975]  and  Gupta  [1976].   The  1975  paper 
uses  the  term  modifier  polynomial   for  A(x).   (In  the  1973  paper  this 
term  is  applied  to  A(x/h).) 

A  very  simple  identification  of  these  polynomial  schemes  with 

Nordsieck  formulas  becomes  apparent  if  the  polynomials  are  represented 

by  vectors  of  scaled  derivatives.   Let  a  =  [p  (t  ),  hp' (t  ),..., 

~n     n  n     n  n 

hqp(q)(t  )/q!]T  and  I   =  [A(0) ,  A'(0),...,  A(q)  (0) /q ! ]T.   Then  from 
n    n  ~ 

(3.1)  we  have 
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p(j)(t  )  =  P(j1)(tn)  +  h"jA(j)(0)6  , 
n    n     n-±  n  n 


from  which  it  easily  follows  that 


a  =  Pa  ,  +  16 
~n    ~n-l   ~  n 

-IT  T 

with  6   chosen  so  that  h  e,a  =  f(t  ,  e.a  ). 
n  ~l~n      n  ~0~n 

For  the  k-th  order  backward  differentiation  formula 

k 

Z  a.y   ,  =  hy' 
j-0  3  n"J 

the  modifier  polynomial  A(x)  must  satisfy 

A(-j)  =  0  ,    j  =  l(l)k-l  , 
A(0)  =  1,  A'(0)  =  aQ  , 


whence 


k-1 
A(-k)  =  (A'(0)  -  Z  a.A(-j)}/a, 


j=0  J 


=  0  . 


Therefore 


A(x)  =  (X^k) 

3=0  k!   ^+1 

where  the  (square)  brackets  denote  Stirling  numbers  of  the  first  kind 

(see,  for  example  Knuth  [1968,  p.  66]).   The  fact  that  A(x)  vanishes 

for  x  =  -1 ,  -2 , . . . ,  -k  means  that  the  values  y  , ,    y     „,...,y  ,  are 

n-1   n-z       n-k. 

"remembered"  after  advancing  from  t   ,  to  t  . 

n-1     n 

For  the  (k+l)-st  order  Adams-Moulton  formula 
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y  -y     =  h  I      B.y'  . 
n   n_1      j=o   J  n~J 

the  modifier  polynomial  of  the  (k+2)-value  form  must  satisfy 

A'(-j)  =  0  ,    j  =  2(l)k-l  , 

A(0)  =  3Q  ,     A»(0)  =  1  , 


whence 


A(-l)  =  A'(-l)  =  0  , 


k-1 

A'(-k)  =  {A(0)  -  A(-l)  -  Z  3.  A,(-j)}/31 


J-0  J 


=  0 


Therefore 


and  so 


A»(x)  =  (X+k) 


A(x)  -  /  (7tk)dy 
-1   k 


=  -/(y>y  +  /  (y>y  . 

0  0 

In  terms  of  Stirling  numbers 

/  C k  )dy  =   Z  TTTT  [  ,  ]xJ  . 

o   k     j=i  Jk'   J 

These  are  the  formulas  used  by  the  nonstiff  options  of  the  codes  DIFSUB 
of  Gear  [1971],  GEAR  Rev.  3  of  Hindmarsh  [1974],  and  EPISODE  of  Byrne 
and  Hindmarsh  [1975]  when  the  stepsize  does  not  vary. 

The  (k+l)-st  order  Adams-Moulton  formula  could  also  be  written 
as  a  (k+l)-value  Nordsieck  formula.   The  modifier  polynomial  for  such  a 
formula  can  be  determined  by  applying  the  theorem  that  follows. 
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THEOREM  3.1.  Let   A(x)  be  the  modifier  polynomial  of  the   (q+1)- 
value  form  of  a  linear  k-step  formula  of  order  at  least   q  where 
q  >  k+1  and   p(£)  and   a(£)  have  no  common  factors.      Then   A(x)  :=  A(x)  -  A(x-l) 

is  the  modifier  polynomial  for  the  q-value  form  of  the  multistep 
formula. 

Proof.      For  m  <^  k 

L  mA(-m)  =  L  A(-m)  -  Z     (-a.A(-i)  +  0.A'(-i)} 

1=0    1        X 

=  0  , 
and  for  m  =  k+1 

30        1   k 
A(0)  =  -f-  A'  (0)  +  -f-  Z  {-a.A(-i)  +  a.A'(-i)} 

ao       ao  i=l    1       X 
=  V 

Hence 

L:JA(-m)  =  0  ,    j  =  0(l)k-m  , 

A(0)  =  3Q,       A'(0)  =  aQ  , 

A  (-j)  =  A'(-j)  =0,    j  =  l(l)m-l  , 
from  which  the  theorem  follows. D 

Therefore  the  (k+l)-value  form  of  the  (k+l)-st  order  Adams- 
Moulton  formula  has 

A(x)  =  /  (y+k)dy  -  f    (^k)dy 
-1   R       -1    k 

=  /  (yf)dy  +  /  {(^k)  -  (y+^1)}dy 
-1   K       0 


■  -  "'  C^ +  {  ?Z>y  ■ 
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These  are  the  original  formulas  of  Nordsieck. 

Let  us  consider  now  the  case  m=0.   The  values  <5  ,  y  ,  and  y' 

n   n       n 

are  determined  as  before.   However,  for  j  =  0(l)k-l 

I/}d  (t  )  =  -a,  .(y  -y  _)  +  hB,  .  (y'-y'  n) 
h  n  n      k-j   n  n,0      k-j   n  n,0 

=  (-a,  .Bn  +  8,  .cx.)6   . 
k-j  0    k-j  0  n 

Therefore  the  solution  is  advanced  one  step  as  given  by  (3.1)  with  A(x) 

determined  by 

L^A(O)  =-ak_.  3Q  +  3k_.aQ,    j  =  0(l)k-l. 

However,  p  (t)  does  not  necessarily  interpolate  y  and  y1  nor  does 
n  n      n 

it  generally  satisfy  the  differential  equation  at  t  =  t  .   Therefore, 

in  the  scaled  derivative  form  this  scheme  is  a  generalized  Nordsieck 

formula  in  the  sense  that  6   is  determined  by  a  condition  other  than  the 

n 

satisfaction  of  the  differential  equation  at  t  =  t  .   Such  formulas  are 

n 

potentially  useful  because  of  their  minimum  storage  property.   Nordsieck 

[1962,  p.  27]  considers  the  possibility  of  having  Z     =  $  but  JL  f   a 

so  that  p  (t)  interpolates  y  but  not  y' ;  however,  the  results  of  his 
rn  n  n 

experiments  were  not  promising.   Also,  Wallace  and  Gupta  [1973]  mention 
that  "Other  choices  of  6   are  rationally  possible,  and  we  hope  to  explore 
some  other  choices  later."  Finally,  it  is  worth  noting  that  Theorem  3.1 
extends  to  the  case  q  =  k. 
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4.   The  Correspondence  between  Multistep  and  Nordsieck  Formulas.   In  the 
preceding  section  it  was  shown  how  to  construct  the  modifier  polynomial 
A(x)  of  a  Nordsieck  formula  from  the  polynomials  p(£)  and  o(E,)    of  a  linear 
multistep  formula  provided  that  p(£)  and  a(£)  have  no  common  divisors. 
In  this  section  we  show  how  to  obtain  p(£)  and  o(E,)    from  A(x),  and  we 
establish  a  one-to-one  correspondence  between  (i)  the  class  of  all  (q+1)- 
value  linear  Nordsieck  formulas  and  (ii)  the  class  of  all  linear  multistep 
formulas  of  formal  order  at  least  q  and  of  stepnumber  at  most  q  (including 
those  for  which  p(£)  and  o(E,)   have  common  factors). 

For  each  linear  Nordsieck  formula  we  define  a  corresponding 
linear  multistep  formula  (p,  a)  by 


(4.1) 


p(0    :=  det(Cl-P)e],(?I-P)  h, 


5(0  :=  det(£I-P)eJ(£I-P)  1%. 


T      -1 
Applying  Cramer's  rule  to  e.(^I-P)   £  for  j  =  1  and  j  =  0  yields 


P(£)  =  det 


5-1 

o 

0 


-1 


-2 


-1 


and 
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5(0  -  det 


% 

-l 

-1 

-l 

*1 

?-i 

-2 

-2 

h 

0 

5-1 

-(q) 

• 
• 
• 

• 
• 
• 

£ 
.  q 

0 

0 

5-1 

Hence  5(C)  is  a  polynomial  of  degree  q  or  less,  and  p(C)  is  of  degree  q 
since  £  ^  0.   Let 


p(5)  =:   Z  a.?q  1  and  5(C)  =:   Z  3.Cq  1 


i=0 


i=0 


It  might  happen  that  a  =3  =0,  and  hence  express 

q   q 


(4.2) 


p(5)  =:  Cm  1P(?)  and  5(C)  =:  C™  Vo 


2    2 
where  a,  +  3,  >  0  and  m-1  =  q-k. 
k    k 

The  following  theorem  shows  that  the  p(C)  and  a(C)  corresponding 
to  A(x)  can  be  used  to  reconstruct  A(x)  by  the  process  of  §3  if  it  is 
applicable. 

THEOREM  4.1.  With   P(C)  and   a  (5)  derived  from   A(x)  by    (4.1)  and 
(4.2)  we  have 


L^A(-m)  =0,j=  0(l)k-l  , 


£Q  =  A(0)  =  3Q  ,  %1   =  A'(0)  =  aQ 


A(-j)  =  A'(-j)  =  0  ,  j  =  l(l)m-l  , 


A(-m)  =  (-Dq3k  ,  A'(-m)  =  (-1)\  . 
Proof.      The  relations  £  =  3n  and  £  =  a  follow  from  the 


expressions  for  p(C)  and  a(C)  as  determinants.   From  (4.1)  we  have 
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P(U    =  -(?-l)q+1e^(I-?P   X)    h   X£ 

oo 

=  -(?-Dq+1  z  Jr^p-j-^ 

j=o  ~x 

CXD 

=  -(?-l)q+1     Z     A'(-j-l)?j    . 
J-0 

Clearly,   A'(-j)   =  0,    j    =   l(l)m-l,   ak  =    (-l)qA'(-m),    and 

00 

p(0    =  -(?-l)q+1      Z     A'(-j-m)CJ    . 
j=0 

Similarly  A(-j)  =  0,  j  =  l(l)m-l,  3R  =  (-l)qA(-m),  and 

oo 

a(C)  =-(C-Dq+1  2  A(-j-m)?j  . 
3=0 

From  these  expressions  for  p(£)  and  a(£)  we  get 

OO  00 

-p(0   Z  A(-j-m)?j  +  o(0      Z  A'(-j-m)?j  =  0  . 
3=0  j-0 

Equating  coefficients  of  the  powers  of  £  completes  the  proof .D 

The  next  theorem  shows  that  the  multistep  formula  (p,  5) 

corresponding  to  A(x)  is  of  order  at  least  q,  and  therefore  it  can  be  used 

to  reconstruct  A(x)  by  the  process  of  §3  if  p(£)  and  a(£)  have  no  common 

factors. 

THEOREM  4.2.  The   linear  multistep  formula  corresponding  to 

a   (q+1) -value   linear  Nordsieck  formula  has  formal  order  at  least   q. 

T 
Moreover j   it  is  of  order  at  least   q+1  if  and  only  if  b  I   =  0  where  the 

components  of  b  are  the  Bernoulli  numbers  defined  by 

00  b.zJ 
eZ-l    j-0  J' 
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Proof.      According  to  the  definition  of  formal  order  and 
equation  (4.1)  it  must  be  shown  that 

zq+1e[((l+z)I-P)_1£  =  log(l+z)zq+1eQ((l+z)I-P)"1£  +  0(zq+1)  , 

and  so  it  is  enough  to  show  that 

e^(I-z_1F)-1  =  log(l-fz)eJ(I-z"1F)~1  +  0(z) 

where  F  :=  P-I.   Because  F  "  =  0,  it  is  not  difficult  to  verify  that 

log(l+z)(I-z_1F)-1  =  log(I+F)(I-z~1F)~1  +  0(z)  , 

T  T 

and  so  it  suffices  to  show  that  e_  log(I+F)  =  e..  .   This  holds  because 

I+F  =  P  =  exp(Q)  where 


Q  = 


o 


o 


o 


To  prove  the  second  assertion,  retrace  the  first  two  steps  of  the  proof 
to  get  the  following  condition  for  formal  order  of  at  least  q+1: 

eJ{l-|F+...+  ^r(-F)lH  =  0  . 

The  matrix  in  braces  is 

log  (i+F)  =   Q    =  2  j_  nJ 

F      '  exp(Q)-l    ±10  V-  ' 


from  which  the  second  part  of  the  theorem  follows. □ 

Remark   1.   The  matrices  P,  F,  and  Q  obviously  represent  the 
linear  operators  shift,  forward  difference,  and  derivative,  respectively, 
for  polynomials  of  degree  q  or  less.   Defining  B  =  I-P    to  be  the  matrix 
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representing  the  backward  difference  operator  V  it  is  not  difficult  to 

T 
show  that  the  condition  b  Z   =  0  is  equivalent  to 

q   1   1 

Z     ~r   VJA(-1)  =  0 

j=0  JT± 

as  well  as 

q-1      , 
A(-l)  -  Z     y^.1VJA,(-1)  =  0 
j=0  3+i 

where  the  y*.,    are  the  coefficients  for  the  backward  difference  form 

of  the  Adams-Moulton  formulas.   These  conditions  on  the  polynomial  A(x) 

arise  in  another  situation  in  Henrici  [1962,  p.  342]. 

T 
Remark   2.   The  condition  b  I   is  obtained  by  Wallace  and 

Gupta  [1973]  although  the  coefficients  b.  are  not  identified  as  the 

Bernoulli  numbers.   Instead  a  recursive  definition  is  given  for  the 

b . ,  which  should  read 

j    1+1 
b  =  1,  Z      (37)b.  =  0  . 

i=0 

Remark   3.   It  is  shown  by  Gear  [1966]  that  a  stable  linear 

Nordsieck  method  is  convergent  of  order  q  in  all  components  of  a  ,  and 

~n 

it  is  shown  by  Skeel  and  Jackson  [1977]  that  it  is  convergent  of 

T 
order  q+1  in  all  components  if  and  only  if  b  Z  =   0. 

It  has  thus  been  shown  that  to  each  linear  Nordsieck  formula, 

which  is  uniquely  specified  in  terms  of  £,  there  corresponds  a  linear 

multistep  formula  of  formal  order  at  least  q  and  of  stepnumber  at  most  q, 

We  now  establish  a  correspondence  in  the  opposite  direction.   We  have 

from  (4.1)  that 
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5(l+z)  -  zq+1eJ(zI-F)  1£ 

=  zqeT     I      Z"jFj£ 
J-0 

where  F  :=  P-I  is  the  matrix  of  the  forward  difference  operator.   Equating 

coefficients  of  powers  of  z  gives 


2^  =  A^A(O)  . 


and  hence  A(x)  is  determined  from  its  corresponding  5(0  by  Newton's 
forward  difference  formula 

q  a(j)Q)   x 
(4.3)  A(x)  =  I  Q  ./   (  X.)  . 

3=0   J!     q^ 
Since  the  linear  multistep  formula  is  uniquely  specified  by  5(?) ,  the 
one-to-one  correspondence  is  established.   The  inverse  transformation 
is  conveniently  expressed  as 

q 

a(U  =  Z  AJA(0)(C-Dq  J 
j=0 

and 

q-1   . 
p(?)  =   E  AJA'(0)(^-l)q  J  . 
j-0 

Note  from  (4.3)  that  the  popular  normalization  a(l)  =  1  corresponds  to 
£  =  1/q!. 

q 

Remark.      Theorem  3.1  can  be  generalized  to  any  modifier 
polynomial  A(x)  for  which  m  >_  2.   Express  d(£)  =:  ^d(^).   Then 

a(j)d)  =  a(j)d)  +  ja(j-1}(i)  , 

from  which  it  can  be  shown  that 
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q_1  S(j)(l)   x+1 
A(x)  =  I  .,U   (X L) 

J-0  j!  q'J 


and 


q-i  a(j)d)       x 

A(x)  :=  A(x)  -  A(x-l)  =  I     . )L)    (     *  ,)  . 

j-o       j!       q_1-J 

Clearly  A(x)  is  the  modifier  polynomial  for  the  q-value  formula, 
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5.   Equivalence  of  Linear  Nordsieck  Formulas  to  Linear  Multistep  Formulas, 
Recall  that  a  linear  Nordsieck  formula  determines  successive  values  by 
the  system  of  difference  equations 

(5.1)  a  =  Pa   .  +  £6 

~n    ~n-l   ~  n 

where  6   is  chosen  so  that  y'  =  f(t  ,  y  ).   We  do  not  consider  the  more 
n  n      n   n 

general  formula  a  =  Aa   ,  +  Z6     because  formulas  with  A  4-   P  are  of 
~n    ~n-l   ~  n 

dubious  value,  and  in  any  case,  most  of  the  results  for  A  =  P  generalize 

if  minor  restrictions  are  imposed  on  A. 

For  theoretical  purposes  a  rewriting  of  (5.1)  is  often  useful. 

Premultiplying  (5.1)  by  e?  yields  6  =  £~1(hf  -e^Pa   ,).   Let  I    :=  C^l 

J   ~1         n    1    n~l  ~n-l        ~     1  ~ 

~  x 
and  S  :=  (I-£e  )P,  and  we  get 

(5.2)  a  =  Sa  ,  +  h£f 

~n    ~n-l    ~  n 

T 
where  f   is  chosen  so  that  f  =  f(t  ,  e.Sa  -  +  h£_f  ).   Expressions  for 
n  n      n  ~1  ~n-l     On 

p(£)  and  5(5)  in  terms  of  S  are  given  by  the  following  theorem,  whose 
proof  uses  an  idea  from  Osborne  [1966,  equation  (4.3)]. 

THEOREM  5.1.  The  polynomials   p(£)  and   5(£)  defined  by    (4.1) 


satisfy 


and 


p(£)  =  det(?I-S)e^(?I-S)_1£ 


5(0  =  det(€l-S)eJ(£I-S)  ^  . 


Proof.      We  have 


£l-S  =  £I-P  +  Ze^P 


=  (?I-P)(I  +  (Cl-P)  ^P)  , 
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and  so 

det(£l-S)'(5l-S)-1£  =  det(Cl-P)  det  (I+C^I-P)"1^?) 

•(I+(?I-P)"1£e^P)"1(Cl-P)"1I  . 

This  can  be  simplified  by  means  of  the  identity 

T       T  -1 
det(I+xy  )'(I+xy  )   x  =  x  , 

which  follows  from  Cramer's  rule.   Thus  we  get 

det(£l-S)'(£l-S)-1£  =  det(5l-P)*(?I-P)_1^  , 
from  which  the  lemma  immediately  follows. □ 

Note  that  e^S  =  0T;  and  so  p(£)  -  det (^I-S)^"1^  ,  which 
gives  the  characteristic  polynomial  of  S  as 

det(£I-S)  =  £~Vp(0  • 

Thus  the  strict  root  condition  is  satisfied  by  the  linear  Nordsieck  formula 
(cf.  Skeel  and  Jackson  [1977])  if  and  only  if  it  is  satisfied  by  the 
corresponding  linear  multistep  formula. 

The  theorem  that  follows  shows  that  in  the  case  of  all   Nordsieck 
formulas  the  zero-th  and  first  components  of  the  vectors  a  satisfy  the 
difference  equation  of  the  corresponding  linear  multistep  formula  (p,  o) . 
Hence  all  the  limitations  on  the  accuracy  of  multistep  formulas  (Dahlquist 
[1956],  [1963])  apply  also  to  Nordsieck  formulas. 

THEOREM  5.2.  The  values   y  ,  y' ,  n  -  0(1)N,  computed  from   (5.1), 

satisfy  the   linear  multistep  formula 

p(E)y    =  ha(E)y' 

n-k         n-k 

defined  by    (4.1)  and   (4.2)  for   n  =  k(l)N. 

Proof.      From  (5.2)  we  have 
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,  •  k-1   . 

Sk_:ia  ,  =  a   .  -   Z   S1  Jilhf   . 
~n-k   ~n-i    .  .     ~   n-i 
i=J 

and 

k  k-1    k-1   ._.- 

p(S)a  ,=  Z  a. a   .-  Z  a.   Z  S  J£hf 

~n-k    .  _  j~n-j    s  n  J  •  j     ~   n-1 
j=0  J    J    j=0    i=j 

T 
Premultiplying  by  e_  and  rearranging  gives 


k  k-1     i         _ 

Z  a.y   .  =  Z  el{  Z  a.S  J  Hhf   .  +  e„p(S)a  .  . 
3=0  *  °-J    1=0  ~°  j-0  3     ~  n_1   ~°    ~n_k 

It  needs  to  be  shown  first  that 

q      i      •  • 
(5.3)  a(0  =  Z  e^{  Z  a.S1  JH£q  X  . 

i=0  ~U  j-0  J 

We  have  that 

=  pCOeJCl-r^)"1! 

=  Z  a.E,q   3   Z  K     enSXi   , 
j=0  3     i=0    ~   ~ 

which  leads  to  (5.3).   It  remains  to  be  shown  that 

Because  S  satisfies  its  characteristic  polynomial,  e  p(S)S  =  0  .   Also, 
because  S  is  of  rank  q,  it  has  only  one  linearly  independent  left  eigenvector 

associated  with  the  eigenvalue  0;  and  so  enp(S)S   '  is  a  multiple  of  e..  . 

T     m-1    ~  T 
This  implies  e  p(S)S    (I-let)    =  0,  and  so  by  (5.3) 

T  ,_,,._m-l    n  T 

?0P(S)S    =  ek+m-l!l  • 
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If  m  >  1,  then  e  p(S)S    =  0  ,  and  this  argument  can  be  repeated  until 

the  assertion  is  established. □ 

Remark   1.   We  have  y1  =  f   except  possibly  for  n  =  0. 

n    n 

Remark   2.   This  theorem  improves  the  result  of  Descloux  [1963] 
and  Osborne  [1966]  in  two  respects:   first,  p(?)  and  cr(£)  are  given  in 
closed  form,  and  second,  the  result  is  shown  for  n  >.  k  rather  than  n  >^  q. 

The  next  theorem  proves  that  Nordsieck  and  multistep  methods 

are  equivalent  in  the  sense  of  Gear,  for  it  is  shown  that  there  is  a 

nonsingular  matrix  that  relates  the  scaled  derivative  approximations  a 

~n 

to  certain  linear  combinations  of  computed  y  and  y'  values.   As  a 
consequence,  in  practically  all  cases  linear  Nordsieck  methods  cannot 
be  regarded  as  generalizations  of  linear  multistep  formulas  (cf.  Gear 
[1971a,  pp.  102,  136]). 

THEOREM  5.3.  Assume  the  polynomials   p(£)  and  o(0   corresponding 
to  a  given  (q+l)-value  Nordsieck  formula  have  no  common  factors.      Then 
there  exists  a  unique  nonsingular   (q+l)x(q+l)  matrix   T  depending  only  on 
H  such  that 

a  =  T_1y 
n      yn 

for   n  =  k-l*(l)N  where 

r  0         k-m-1*  ,  ,      ,  f     -.T 

y   :=  [s    ,...,  s      ,  y   *,...,  y     ,  hy  , .  .  .  ,  hy     J 
in      n-m       n-m      n-0        n-m+1    n        n-m+1 

Proof.      By  Theorem  4.2  the  formula   (p,  a)  is  of  order  at  least 

q,  and  hence  by  the  corollary  to  Theorem  2.1  there  exists  a  unique 

nonsingular  matrix  T  such  that 


31 


P(tn) 


hp'(tn) 


hS(q)(t  )/q! 


=  T 


-1 


LTp(t    ) 
h   n-m 


L"p(t    ) 
h   n-m 


"^Wi' 


for  any  polynomial  of  degree  q  or  less.   Define  values  y.  and  y!  for 
j  =  l*-k(l)-l  such  that 

The  process  of  §3  may  be  applied  to  construct  a  (q+l)-value  linear 

Nordsieck  formula  which  would  compute  vectors  T  y  ,  n  =  1(1)N.   The 

results  of  §4  imply  that  this  formula  would  be  identical  to  the  given 

Nordsieck  formula,  and  hence  a  =  T  y  .   The  uniqueness  of  T  follows 

~n      in 

from  the  fact  that  the  method  is  exact  for  all  polynomials  of  degree 

at  most  q  if  the  starting  values  are  exact. □ 

Remark.      For  implicit  formulas  Wallace  and  Gupta  [1973]  give 

an  informal  argument  suggesting  that  the  quantities  a    ,  and  6   . , 

~n-q       n-j 

j  =  l(l)q-l,  can  be  expressed  as  linear  combinations  of  y   .,  f   ., 

n-j   n-j' 

■m—  1 

j  =  l(l)q.   Their  conjecture  is  correct  as  long  as  E,        p(K)    and 

E,        O(C)  have  no  common  factors.   In  that  case  it  is  an  immediate 

consequence  that  the  components  of  a    are  linear  combinations  of 

y   .,  f   .,  j  =  l(l)q. 
n-j    n-j 

Theorem  5.3  does  not  extend  to  Nordsieck  formulas  for  which  the 
corresponding  p  and  G  polynomials  have  common  factors. 

THEOREM  4.4.  Assume  that   p(£)  and   a(£)  have  a  common  factor 
and  that  the  formula   (p,  a)   is  of  order  at   least   q.  Then   a  cannot  be 


expressed  only  in  terms  of  y   .,  hy'  .,  j  =  0(l)n. 
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Proof.      We  begin  by  showing  that  there  exists  an  eigenvector 

T 
v  associated  with  a  nonzero  eigenvalue  £,  of  S  such  that  e~v  =  0.   From 

(4.3)  it  follows  that  p(£.)  and  a(£)  have  the  common  factor  1  if  and  only 

if  £  =0.   First,  suppose  that  £  =   0.   Then  £  =  1  is  a  common  root  and 

q  q 

T 
so  the  formula  must  be  of  formal  order  at  least  k+1.   Hence  b  £  =  0 

where 

T    T  log(P) 
H    fo   P-I   ' 

and  so 

rp     rp  rp  rp  rp  rp  rp 

(b  +e|)S  =  b  P  =  b  +  ej  logP  -  b  +  e£  . 

T     T 
Also  £  =  0  implies  e  S  =  e  .   Therefore  the  null  space  of  S-I  is  at  least 
q  ~q     ~q 

of  dimension  two,  and  so  we  can  choose  v  ^  0  such  that  Sv  =  v  and 

T 
e.v  =  0.   Second,  suppose  that  £  ^  0.   Let 
~0~  q 

v(K)    =   det(£l-S)-(a-S)_1£  , 

which  is  a  vector  of  polynomials  in  £.   Let  £,  be  a  common  root  of  p(0 

and  a(£).   Then  v(£_)  ^  0,  since  e  v(£)  =  (£-l)q£  .   Also 

(J     ~        ~q~  q 

(C0I-S)y(£Q)  =  det(£QI-S).£  =  0  , 

T 
and  env(E.  )  =  a(£.  )  =  0.   Therefore  in  either  case  there  exists  E,     £   0 

T  T 

and  v  4-   0  such  that  S  =  £_v  and  e~v  =  0.   Moreover,  e,v  =  0,  because 
~  v    0~     ~0~  ~1~ 

T 
e..  is  a  left  eigenvector  for  the  eigenvalue  0.   This  means  that  if  a_ 

were  changed  to  aA  +  v,  the  values  of  y  and  hy'  would  remain  unchanged 
~U   ~  n       n 

for  all  n,  and  yet  a  would  become  a  +  K„v.      This  proves  the  nonexistence 

~n  ~n    0~ 

of  an  expression  for  a  in  terms  of  the  values  y  and  y1 .D 

n  n      n 

Theorems  5.3  and  5.4  do  not  cover  the  case  where  the  order  is 
less  than  q  due  to  a  common  factor  of  £-1,  but  extending  the  results  to 
such  formulas  does  not  seem  worthwhile. 
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6.   Equivalence  of  Predictor-corrector  Formulas.   The  use  of  a  linear 

Nordsieck  formula  requires  the  solution  for  <5   of  the  equation 

n 

hy'  n   +  JL6  =  hf(t  ,  y  n  +  1.8   ) 
■'njO    In       n  ^n,0    0  n 

T  T 

where  y   _  =  e^Pa   ,  and  hy1  .  =  e.Pa   ,  .   In  practice  this  must  be 
'n,0   ~0  ~n-l       n,0   ~1  ~n-l 

approximately  solved  by  some  iterative  process: 

6  n  :-  0  , 

n,0 

6    ..  :=  6   -[£l-£hJ   ]~1[hy*  _+£_  6   -hf   ]  , 
n,y+l     n,y   1   0  n,y      n,0  1  n,y   n,y 

M  =  0,  1,...,  M(n)-1  , 

where  f    ,1,  and  J    ,  are  defined  as  in  §2  with  y    =  y   _  +  &„6 

n,y         n,y-l  n,y    n,0    0  n,y 

After  determining  8      :=  6  .,,  N  and  adding  the  increment  £6   a  final 

n     n,M(n)  ~  n 

function  evaluation  may  or  may  not  be  performed.   For  a  P(EC)*  Nordsieck 

formula 

a   :=  Pa  ,  +  18 
~n     ~n-l    ~  n 

and  for  a  P(EC)*E  Nordsieck  formula 

a   :=  Pa  .  +  £6  +  e.h(f  -y ' )  , 
~n     ~n-l   ~  n   ~1   n  n 

where  hy'  :=  hy'  .   +  1.8      „,  ,  n . 
n     yn,0    1  n,M(n)-l 

Remark.      The  predictor-corrector  Nordsieck  formulas  of  Gear 

[1971a]  are  introduced  independently  of  linear  ("corrector  only") 

Nordsieck  formulas.   For  this  reason  these  predictor-corrector  formulas 

use  l^l-lnhJ  :=  1. 

1   0  n,y 

T 
For  the  P(EC)*  formulas  we  have  that  hy'  =  e„Pa   n  +  16  , 

n   ~1  ~n-l    1  n 

and  so  the  recurrence  can  also  be  expressed  as 

a  =  Sa   .  +  £hy'  . 
~n    ~n-l   ~  n 

Therefore  the  equivalence  results  apply  to  P(EC)*  formulas  as  much  as 

they  do  to  linear  formulas.   The  underlying  multistep  predictor  formula 
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T 
can  be  obtained  from  y   .  =  e^Pa   ,  by  expressing  a   n  as  linear 

n,0   ~0  ~n-l  ~n-l 

combinations  of  the  components  of  y  _. ;  thus  it  has  stepnumber  of  at  most 

k+0*. 

The  situation  is  quite  different  for  P(EC)*E  formulas. 

To  determine  the  equivalent  multistep  formula,  one  begins  with  the 

recurrence 

a  =  Sa  ,  +  (i-ejhy1  +  e,hf   . 
~n    ~n-i    ~  ~1   n   ~1  n 

Proceeding  as  in  the  proof  of  Theorem  5.2,  one  obtains  a  difference 

equation  involving  y   . ,  j  =  0(l)k,  y1  . ,  j  =  0(l)k-l,  and  y'  ., 

j  =  l(l)k,  which  is  not  a  true  P(EC)*E  multistep  formula. 

For  example,  consider  the  three-value  Nordsieck  formula. 

By  expressing  y  and  y   .  as  functions  of  y1  ,  y'  ,  ,  f   .,  ,  and  a   _ 
n      n-1  n   n-1   n-1      ~n-2 

and  then  eliminating  h  y"  0/2,  one  obtains  the  difference  equation 

n-z 

\   +  (2,2-2)Vl  +  (1-UlK-2 

-   40h"yn  +  [(1^0)hy;-l  +  ^2-£0)hCl]  +  (%+£2"1)hy;-2  > 
which  is  not  a  P(EC)*E  formula  unless  £  =  £  .   For  the  third  order 
Adams-Moulton  formula  this  is 

y„  -  Vi  ■  Ti  K  +  lT2  K-i +  n  hK-i]  +  s  hy;-2  • 
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7.   Applications.   An  important  practical  consequence  of  the  equivalence 
theory  is  that  all  multistep  formulas  are  minimum  storage  formulas. 
This  idea  is  certainly  implicit  in  the  investigations  by  Gupta  and 
Wallace  of  new  multistep  formulas.   It  is  also  the  rationale  for  the 
computer  research  of  Kong  [1977]  for  k-step  formulas  of  order  k  having  the 
smallest  error  coefficient  for  a  given  angle  a  of  A(a)-stability.   Previous 
searches  by  Dill  and  Gear  [1971]  (the  error  coefficient  for  their  formula 
should  be  14.0  rather  than  0.0108)  and  by  Jain  and  Srivastava  [1970] 
concentrated  on  formulas  for  which  most  of  the  trailing  coefficients  of 
a(£)  were  set  to  zero.   The  computer  results  of  Kong  suggested  formulas 
which  lead  to  a  proof  of  the  following  result:   for  any  a  <  tt/2  there 
exists  an  A(a)-stable  k-step  formula  of  order  k.   This  is  also  proved 
by  Grigorieff  and  Schroll  [1977]. 

The  paper  of  Nordsieck  mentions  that  "the  potential  advantage 
of  a  more  elaborate  procedure  in  which  the  matrix  hf   is  numerically 

y 

computed  at  every  step  and  £  is  made  a  chosen  function  of  hf  ,  implying 
a  nonlinear  process  tailored  to  the  subject  differential  equation  system, 
is  an  interesting  topic  for  future  investigation,"  where  we  have 
substituted  our  notation  for  his.   The  idea  is  developed  further  in  a 
report  of  Gear  [1966],  which  proposes  a  formula  for  scalar  implicit 
differential  equations  of  any  order.   For  the  equation  y'  -  f(t,  y)  =  0 
the  formula  becomes 

I    :-  iVhf 

y 


and 


£,l-JlnhJ    :=  l+h2f2 
1   0  n,u        y 
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where  AM  and  BD  refer  to  the  (q+l)-st  Adams-Moulton  and  q-th  order 
backward  differentiation  formulas,  respectively.   (The  subscripts  y 
in  the  equation  on  page  21  of  this  report  should  read  a  . )   This  idea 

q 

is  extended  to  systems  of  differential  equations  by  Skeel  and  Kong  [1977], 
who  suggest 

y 

and 

ItI-IAiJ  :  =  (1-chf  )2 

1   0  n,y         y 

where  y  an<3  c  are  free  parameters.   In  §4  it  was  shown  that  p(£)  and  o(E,) 

are  linear  transformations  of  Z,    and  so  a  linear  combination  of  I 

vectors  corresponds  to  that  same  linear  combination  of  the  corresponding 

linear  multistep  formulas.   Thus,  under  the  assumption  of  constant  hf 

y 

the  "blended"  Nordsieck  formula  is  equivalent  to  the  same  blend  of  linear 
multistep  formulas. 
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